/*
 * 碰撞问题：
 * 	在发生碰撞时，并不改变位移，只是改变速度
 * 万有引力问题：
 * 	重点是质量
 * insert:
 * 	只能插入叶子节点，否则会退出
 * 	对于非叶子节点，必须是在插入叶子节点的过程中顺便插入的
 * delete：
 * 	只能删除叶子节点，否则会退出
 * 	对于非叶子节点，必须是在删除叶子节点的过程中顺便删除的
 * 测试：
 * 	质量——测试通过
 * 	八叉树——测试通过
 * 	insert/delete——测试通过
 * 万有引力的细节性问题：
 * 	八叉树每个节点的重心和半径的更新——测试通过
 **/

#include <stdio.h> 
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <omp.h>

#define L 1000000
#define RMAX 100
#define RMIN 50
#define VMAX 30
#define p 100000000
#define G 6.67*pow(10,-11)

typedef struct NODE
{	
	int i;					//作为孩子节点的编号
	double minx;			//节点x轴范围——最小值
	double maxx;			//节点x轴范围——最大值
	double miny;			//节点y轴范围——最小值
	double maxy;			//节点y轴范围——最大值
	double minz;			//节点z轴范围——最小值
	double maxz;			//节点z轴范围——最大值
	int isPoint;			//是否为叶子节点
	int childN;				//内部节点的孩子节点数量

	int id;					//叶子节点的编号
	double x;				//叶子节点的x方向坐标
	double y;				//叶子节点的y方向坐标
	double z;				//叶子节点的z方向坐标
	double r;				//叶子节点的半径或内部节点的半径
	double vx;				//叶子节点的x方向速度
	double vy;				//叶子节点的y方向速度
	double vz;				//叶子节点的z方向速度
	double ax;				//叶子节点的x方向加速度
	double ay;				//叶子节点的y方向加速度
	double az;				//叶子节点的z方向加速度
	double m;				//叶子节点的质量或内部节点的半径
	
	struct NODE* child[8];	//内部节点的孩子节点指针
	struct NODE* up;		//节点的父亲节点指针
} Node;

Node *point;
Node tree;
int N;

int retID(Node* root,Node* point)
{
	if(point->x < root->minx || point->x > root->maxx
	|| point->y < root->miny || point->y > root->maxy
	|| point->z < root->minz || point->z > root->maxz
	)
	{
		printf("%d:%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",point->id,point->x,point->y,point->z,root->minx,root->miny,root->minz,root->maxx,root->maxy,root->maxz,point->ax,point->ay,point->az);
		exit(1);
	}
	int i,j,k;
	if(point->x < (root->minx + root->maxx)/2)
		i = 0;
	else
		i = 1;
	if(point->y < (root->miny + root->maxy)/2)
		j = 0;
	else
		j = 1;
	if(point->z < (root->minz + root->maxz)/2)
		k = 0;
	else
		k = 1;
	return i+2*j+4*k;
}

void initNode(Node* root,Node* node,int id)
{
	int i;
	for(i = 0;i < 8;i++)
		node->child[i] = NULL;
	node->isPoint = 0;
	node->childN = 0;
	node->m = 0;
	node->up = root;
	node->i = id;
	node->id = -1;
	if(id/4  == 0)
	{
		node->minz = root->minz;
		node->maxz = (root->minz+root->maxz)/2;
	}
	else
	{
		node->minz = (root->minz+root->maxz)/2;
		node->maxz = root->maxz;
	}
	id = id%4;
	if(id/2  == 0)
	{
		node->miny = root->miny;
		node->maxy = (root->miny+root->maxy)/2;
	}
	else
	{
		node->miny = (root->miny+root->maxy)/2;
		node->maxy = root->maxy;
	}
	id = id%2;
	if(id  == 0)
	{
		node->minx = root->minx;
		node->maxx = (root->minx+root->maxx)/2;
	}
	else
	{
		node->minx = (root->minx+root->maxx)/2;
		node->maxx = root->maxx;
	}
	node->r = 0;
}

void initPoint(Node* point,int i)
{
	point->id = i;
	point->i = -1;
	point->isPoint = 1;
	point->childN = -1;
	point->x = rand()%L;
	point->y = rand()%L;
	point->z = rand()%L;
	point->up = NULL;	
		
	point->r = RMIN+rand()%(RMAX-RMIN);
	point->vx = -VMAX+rand()%(2*VMAX);
	point->vy = -VMAX+rand()%(2*VMAX);
	point->vz = -VMAX+rand()%(2*VMAX);
	point->ax = 0;
	point->ay = 0;
	point->az = 0;
	point->m = 4.0/3*3.14*point->r*point->r*point->r*p;	
}

int collision(Node* point1 ,Node* point2)
{	
	double l = sqrt(pow(point1->x-point2->x,2)+pow(point1->y-point2->y,2)+pow(point1->z-point2->z,2));
	if(l < point1->r + point2->r)
	{		
		double v1,v2;
		v1 = (point1->m*point1->vx-point2->m*point1->vx+2*point2->m*point2->vx)/(point1->m+point2->m);
		v2 = (2*point1->m*point1->vx-point1->m*point2->vx+point2->m*point2->vx)/(point1->m+point2->m);
		point1->vx = v1;
		point2->vx = v2;
		v1 = (point1->m*point1->vy-point2->m*point1->vy+2*point2->m*point2->vy)/(point1->m+point2->m);
		v2 = (2*point1->m*point1->vy-point1->m*point2->vy+point2->m*point2->vy)/(point1->m+point2->m);
		point1->vy = v1;
		point2->vy = v2;
		v1 = (point1->m*point1->vz-point2->m*point1->vz+2*point2->m*point2->vz)/(point1->m+point2->m);
		v2 = (2*point1->m*point1->vz-point1->m*point2->vz+point2->m*point2->vz)/(point1->m+point2->m);
		point1->vz = v1;
		point2->vz = v2;
		
		int k = 0;
		do
		{
			if(k == 0)
			{
				point2->x = (point1->r+point2->r)*1.01/l*(point2->x-point1->x)+point1->x;
				point2->y = (point1->r+point2->r)*1.01/l*(point2->y-point1->y)+point1->y;
				point2->z = (point1->r+point2->r)*1.01/l*(point2->z-point1->z)+point1->z;
				if(point2->x < 0)
					point2->x = -point2->x;
				if(point2->y < 0)
					point2->y = -point2->y;
				if(point2->z < 0)
					point2->z = -point2->z;
				if(point2->x > L)
					point2->x = 2*L-point2->x;
				if(point2->y > L)
					point2->y = 2*L-point2->y;
				if(point2->z > L)
					point2->z = 2*L-point2->z;
			}
			else
			{
				point1->x = (point1->r+point2->r)*1.01/l*(point1->x-point2->x)+point2->x;
				point1->y = (point1->r+point2->r)*1.01/l*(point1->y-point2->y)+point2->y;
				point1->z = (point1->r+point2->r)*1.01/l*(point1->z-point2->z)+point2->z;
				if(point1->x < 0)
					point1->x = -point->x;
				if(point1->y < 0)
					point1->y = -point->y;
				if(point1->z < 0)
					point1->z = -point->z;
				if(point1->x > L)
					point1->x = 2*L-point->x;
				if(point1->y > L)
					point1->y = 2*L-point->y;
				if(point1->z > L)
					point1->z = 2*L-point->z;
			}
			k++;
			l = sqrt(pow(point1->x-point2->x,2)+pow(point1->y-point2->y,2)+pow(point1->z-point2->z,2));	
		}while(l < point1->r+point2->r);
		return 1;
	}
	return 0;
}

void gravity(Node* root,Node* point)
{
	double l = sqrt(pow(point->x-root->x,2)+pow(point->y-root->y,2)+pow(point->z-root->z,2));
	if(root->isPoint && (root->id == point->id || l <= root->r+point->r))
		return;
	if(root->r/l <= 3.14/180 || root->isPoint)
	{
		double a = root->m/(l*l)*G;
		point->ax += a/l*(root->x-point->x);
		point->ay += a/l*(root->y-point->y);
		point->az += a/l*(root->z-point->z);
	}
	else
	{
		int i;
		for(i = 0;i < 8;i++)
		{
			if(root->child[i] != NULL)
				gravity(root->child[i],point);
		}		
	}
}

void delete(Node* point)
{
	if(!point->isPoint)
		exit(2);
		
	/*处理多余节点*/	
	int t = point->i;
	Node *temp1;
	Node* temp = point->up;
	while(temp != NULL)
	{
		if(temp->isPoint)
			exit(2);		
		temp->childN --;
		temp->child[t] = NULL;
		temp1 = temp->up;
		t = temp->i;
		if(temp != NULL && temp->childN == 0 && temp->i != -11)
		{
			free(temp);
			temp = temp1;
		}
		else
			break;
	}
		
	if(point->up != NULL)
	{
		point->up->child[point->i] = NULL;
		point->up = NULL;
	}
	
	point->i = -1;	
}

void insert(Node* root,Node* point)
{	
	if(root->isPoint)
		exit(2);
	if(!point->isPoint)
		exit(2);
	int i = retID(root,point);
	if(root->child[i] == NULL)
	{
		root->childN ++;
		root->child[i] = point;
		point->up = root;
		point->i = i;
	}
	else if(root->child[i]->isPoint == 0)
		insert(root->child[i],point);
	else
	{	
		Node* node = root->child[i];	//保存叶子			
		if(collision(node,point))
		{			
			delete(node);
			insert(&tree,node);			
			insert(&tree,point);
		}
		else
		{
			node->i = -1;
			node->up = NULL;
			root->child[i] = NULL;			
			Node* new = malloc(sizeof(Node));		//新的节点/非叶子	
			initNode(root,new,i);
			root->child[i] = new;
			new->up = root;					//处理新节点	
			insert(new,node);			
			insert(new,point);
		}
	}
}

void updateTree(Node *root)
{
	if(root->isPoint)
		return;
	int i;
	double MX = 0,MY = 0,MZ = 0;
	root->m = 0;
	for(i = 0;i < 8;i++)
	{
		if(root->child[i] != NULL)
		{
			updateTree(root->child[i]);
			root->m += root->child[i]->m;
			MX += root->child[i]->m*root->child[i]->x;
			MY += root->child[i]->m*root->child[i]->y;
			MZ += root->child[i]->m*root->child[i]->z;
		}
	}
	root->x = MX/root->m;
	root->y = MY/root->m;
	root->z = MZ/root->m;
	root->r = 0;	
	if(root->childN == 1)
	{
		for(i = 0;i < 8;i++)
		{
			if(root->child[i] != NULL)
			{
				root->r = root->child[i]->r;
				break;
			}
		}		
	}
	else
	{
		for(i = 0;i < 8;i++)
		{
			if(root->child[i] != NULL)
			{
				double l = sqrt(pow(root->child[i]->x-root->x,2)+pow(root->child[i]->y-root->y,2)+pow(root->child[i]->z-root->z,2));
				if(l > root->r)
					root->r = l;
			}
		}
	}
	int k = 0;
	for(i = 0;i < 8;i++)
	{
		if(root->child[i] != NULL)
			k++;
	}
	if(k == 0)
	{
		free(root);
		root = NULL;
	}
}

void move(Node* root,Node* point)
{
	point->x += point->vx;
	point->y += point->vy;
	point->z += point->vz;
	
	if(point->x < 0)
	{
		point->x = -point->x;
		point->vx = -point->vx;
	}
	if(point->y < 0)
	{
		point->y = -point->y;
		point->vy = -point->vy;				
	}
	if(point->z < 0)
	{
		point->z = -point->z;
		point->vz = -point->vz;
	}
	if(point->x > L)
	{
		point->x = 2*L-point->x;
		point->vx = -point->vx;		
	}
	if(point->y > L)
	{
		point->y = 2*L-point->y;
		point->vy = -point->vy;
	}
	if(point->z > L)
	{
		point->z = 2*L-point->z;
		point->vz = -point->vz;
	}
	
	delete(point);
	insert(root,point);
}

void run()
{
	int i;
	for(i = 0;i < N;i++)
		move(&tree,&point[i]);
	updateTree(&tree);
	for(i = 0;i < N;i++)
	{
		point[i].ax	= 0;
		point[i].ay	= 0;
		point[i].az	= 0;
		gravity(&tree,&point[i]);
	}
	for(i = 0;i < N;i++)
	{
		point[i].vx += point[i].ax;
		point[i].vy += point[i].ay;
		point[i].vz += point[i].az;
	}
}

void init (void)
{
	tree.isPoint = 0;
	tree.childN = 0;
	tree.minx = 0;
	tree.maxx = L;
	tree.miny = 0;
	tree.maxy = L;	
	tree.minz = 0;
	tree.maxz = L;
	tree.up = NULL;	
	tree.i = -11;	//表示root
	tree.id = -1;
	tree.m = 0;
	int i;
	for(i = 0;i < 8;i++)
		tree.child[i] = NULL;
	srand(time(NULL));
	for(i = 0;i < N;i++)
	{
		initPoint(&point[i],i);
		insert(&tree,&point[i]);
	}
}

int main(int argc, char *argv[])
{
	int TIME;
	if(argc > 2)
	{
		sscanf(argv[1],"%d\n",&N);
		sscanf(argv[2],"%d\n",&TIME);
	}
	else
	{
		N = 10000;
		TIME = 10;
	}
	point = malloc(sizeof(Node)*N);
	init();
	double Start,End;
	Start = omp_get_wtime();
	double Time = 0;
	int Round = 0;
	do
	{
		run();
		End = omp_get_wtime();
		Time = End - Start;
		Round++;
	}while(Time < TIME);
	printf("the time is : %lf\n",Time/Round);
	return 0;
}
